/*
 * Cnanowire - computes current-voltage characteristic of bare cylindrical nanowire
 * Copyright (C) 2010  Oka Kurniawan
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public License
 * as published by the Free Software Foundation; either version 2
 * of the License, or (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 *
*/
#include <iostream>
#include <itpp/itbase.h>
#include <itpp/base/timing.h>
#include <omp.h>
#include "deviceinfo.h"
#include "globalconst.h"
#include <boost/math/complex/acos.hpp>
#include "nfile.h"

using namespace itpp;
std::complex<double> imath(0,1);

mat buildHamiltonian(double t0,int NG)
{
	mat retM(NG,NG);
	retM.zeros();
	retM(0,0)=2*t0;
  	retM(0,1)=-t0;
  	retM(NG-1,NG-1)=2*t0;
  	retM(NG-1,NG-2)=-t0;
#pragma omp parallel shared(retM)
  	{
#pragma omp for schedule(dynamic) nowait
  	for (int irow=1; irow<NG-1; irow++)
    {
		for(int jcol=1; jcol<NG-1; jcol++)
		{
	  		if(irow==jcol)
	    	{
	      		retM(irow,jcol)=2*t0;
	      		retM(irow,jcol-1)=-t0;
	      		retM(irow,jcol+1)=-t0;
	    	}
		}
    }
  	} /* end of parallel section */
  	return retM;
}

int main(int argc, char** argv)
{
	CPU_Timer mytime;
	mytime.start();
  	vec zn="2.4048 3.8317 5.1356 5.5201 6.3802 7.0156 7.5883 8.4172";
  	vec kn=zn/DeviceInfo::RNW;
  	vec En=pow(kn,2.0)*GlobalConst::M_HBAR*GlobalConst::M_HBAR/(2*DeviceInfo::MR*GlobalConst::QC);
  	vec xgrid=linspace(0.0,DeviceInfo::LDEV,DeviceInfo::NGRID);
  	double dx=xgrid(1)-xgrid(0);
  	vec Vgrid=linspace(0.0,DeviceInfo::VBIAS,DeviceInfo::NV);
  	double t0=GlobalConst::M_HBAR*GlobalConst::M_HBAR/(2*DeviceInfo::MZ*dx*dx*GlobalConst::QC);
  	mat H0=buildHamiltonian(t0,DeviceInfo::NGRID);
  	vec egrid=linspace(-0.8,1.8,DeviceInfo::NE);
  	double dE=egrid(1)-egrid(0);
  	vec fermi1=zeros(DeviceInfo::NE);
  	vec fermi2=zeros(DeviceInfo::NE);
  	vec gT=zeros(DeviceInfo::NE);
  	vec cur=zeros(DeviceInfo::NV);
  	int vloop, mloop, eloop;
  	for (vloop=0; vloop<DeviceInfo::NV; vloop++)
    {
    	double Vnow=Vgrid(vloop);
      	double mu1=DeviceInfo::EF;
      	double mu2=mu1-Vnow;
      	vec pot=linspace(0.0,-Vnow,DeviceInfo::NGRID);
      	fermi1=1.0/(1.0+exp((egrid-mu1)/DeviceInfo::KBT));
      	fermi2=1.0/(1.0+exp((egrid-mu2)/DeviceInfo::KBT));
      	gT.zeros();
      	cur(vloop)=0.0;
      	for (mloop=0; mloop<DeviceInfo::NMODES; mloop++)
		{
      		vec Emsubx=En(mloop)+pot;
	  		mat H=H0+diag(Emsubx);
	  		vec Tm=zeros(DeviceInfo::NE);
#pragma omp parallel default(shared) private(eloop)
	  		{
#pragma omp for schedule(dynamic)
	  		for (eloop=0; eloop<DeviceInfo::NE; eloop++)
	    	{
	      		cmat sigma1(DeviceInfo::NGRID,DeviceInfo::NGRID);
	      		sigma1.zeros();
	      		cmat sigma2(DeviceInfo::NGRID,DeviceInfo::NGRID);
	      		sigma2.zeros();
	      		cmat gamma1(DeviceInfo::NGRID,DeviceInfo::NGRID);
	      		cmat gamma2(DeviceInfo::NGRID,DeviceInfo::NGRID);
	      		cmat G(DeviceInfo::NGRID,DeviceInfo::NGRID);

	      		std::complex<double> h0barleft=2*t0-(egrid(eloop)+imath*DeviceInfo::ETA)+Emsubx(0);
	      		std::complex<double> h0barright=2*t0-(egrid(eloop)+imath*DeviceInfo::ETA)+Emsubx(DeviceInfo::NGRID-1);

	      		cmat Amat(2,2),Bmat(2,2);
	      		Amat(0,0)=h0barleft;
	      		Amat(0,1)=-t0;
	      		Amat(1,0)=1;
	      		Amat(1,1)=0;
	      		Bmat(0,0)=t0;
	      		Bmat(0,1)=0;
	      		Bmat(1,0)=0;
	      		Bmat(1,1)=1;

	      		cvec lambda(2);
	      		cmat evecs(2,2);
	      		eig(Amat,Bmat,lambda,evecs);
	      		vec ikn=imag(-imath*log(lambda));
	      		ivec idkvec=find(ikn>0);
	      		int idk=idkvec(0);
	      		std::complex<double> c0=evecs(0,idk);
	      		std::complex<double> Tleft=c0*lambda(idk)*1/c0;
	      		sigma1(0,0)=-t0*Tleft;


	      		//std::complex<double> ka=boost::math::acos(1.0-(egrid(eloop)
					  //     +imath*DeviceInfo::ETA-Emsubx(0))/(2.0*t0));
	      		//sigma1(0,0)=-t0*exp(imath*ka);
	      		//std::complex<double> buf=-t0*exp(imath*ka);

	      		gamma1=imath*(sigma1-sigma1.hermitian_transpose());

	      		Amat(0,0)=h0barright;
	      		Amat(0,1)=-t0;
	      		Amat(1,0)=1;
	      		Amat(1,1)=0;
	      		Bmat(0,0)=t0;
	      		Bmat(0,1)=0;
	      		Bmat(1,0)=0;
	      		Bmat(1,1)=1;

	      		eig(Amat,Bmat,lambda,evecs);
	      		ikn=imag(-imath*log(lambda));
	      		idkvec=find(ikn<0);
	      		idk=idkvec(0);
	      		c0=evecs(0,idk);
	      		std::complex<double> Tright=c0*1/(lambda(idk))*1/c0;
	      		sigma2(DeviceInfo::NGRID-1,DeviceInfo::NGRID-1)=-t0*Tright;

	      		//ka=boost::math::acos(1.0-(egrid(eloop)
	      			//+imath*DeviceInfo::ETA-Emsubx(DeviceInfo::NGRID-1))/(2.0*t0));
	      		//sigma2(DeviceInfo::NGRID-1,DeviceInfo::NGRID-1)=-t0*exp(imath*ka);
	      		//buf=-t0*exp(imath*ka);

	      		gamma2=imath*(sigma2-sigma2.hermitian_transpose());
	      		G=inv(((egrid(eloop)+imath*DeviceInfo::ETA)*eye(DeviceInfo::NGRID))
		    			-H0-diag(Emsubx)-sigma1-sigma2);
		    	Tm(eloop)=real(trace(gamma1*G*gamma2*G.hermitian_transpose()));
	    	}
	  		}/*end of parallel section */
	  		gT+=Tm;
	  	}
      	cur(vloop)=2*GlobalConst::I_CONST*sum(elem_mult(gT,(fermi1-fermi2)))*dE;
    }
  	NFile fcur("ivcurve.dat");
  	fcur.addData("Vbias",Vgrid,0);
  	fcur.addData("Cur",cur,1);
  	fcur.saveToFile();
  	mytime.stop();
  	std::cout<<"time (s): "<<mytime.get_time()<<"\n";
  	return 0;
}

